import numpy as np


def Method_Iterative(A_shape, b, D, H, Precision):

	D_Iteration, H_Iteration, Constant_Matrix, Initial_Vector = Generate_Matrix_for_Iterate(A_shape, b, D, H)
	#print(Iteration_Matrix)
	#这里的Ny和之前的Ny是同一个值
	Ny = D_Iteration.shape[0]
	Iteration_Matrix_First = np.bmat([D_Iteration, H_Iteration])
	Iteration_Matrix_Middle = np.bmat([H_Iteration, D_Iteration, H_Iteration])
	Iteration_Matrix_Last = np.bmat([H_Iteration, D_Iteration])
	#这里的Nx和之前的Nx是同一个值
	Nx = int(A_shape[0]/Ny)
	X_Real_Time = Initial_Vector
	Error = Precision + 1 
	count = 0
	while Error > Precision:
		X = np.copy(X_Real_Time)
		X_Real_Time[0:Ny] = Calcule_Part_of_X(Iteration_Matrix_First, Constant_Matrix[:Ny], X_Real_Time[:Ny*2])
		for i in range(1,Nx-1):
			X_Real_Time[Ny*i:Ny*(i+1)] = Calcule_Part_of_X(Iteration_Matrix_Middle, Constant_Matrix[Ny*i:Ny*(i+1)], X_Real_Time[Ny*(i-1):Ny*(i+2)])
		X_Real_Time[Ny*(Nx-1):] = Calcule_Part_of_X(Iteration_Matrix_Last, Constant_Matrix[Ny*(Nx-1):], X_Real_Time[Ny*(Nx-2):])
		Error = np.max(np.abs(X_Real_Time - X))
		count = count + 1
	return X, Error, count

def Calcule_Part_of_X(Iteration_Matrix, Constant_Matrix, Initial_Vector):
	X = Iteration_Matrix*Initial_Vector + Constant_Matrix
	return X

def Generate_Matrix_for_Iterate(A_shape, b, D, H):
	"""
	A = M - N 
	N = M - A
	x = M^{-1}*N*x + M^{-1}*b
	x = Iteration_Matrix*x + Constant_Matrix
	"""
	"""
	把 Iteration_Matrix像对A一样分块， 这样可以分块的运算X
	"""
	Dimention = D.shape[0]
	M_inv = 1/D.item(0,0)
	D_Iteration = np.matrix(np.eye(Dimention,Dimention)) - M_inv*D
	H_Iteration = -M_inv*H
	Initial_Vector = np.matrix(np.zeros([A_shape[0],1]))
	Constant_Matrix = M_inv*b
	#D_Iteration, H_Iteration是(Ny,Ny)维，Constant_Matrix, Initial_Vector是(Nx*Ny,1)维
	return D_Iteration, H_Iteration, Constant_Matrix, Initial_Vector


def Generate_Matrix_A_b(Nx, Ny, Xmax, Xmin, Ymax, Ymin):
	"""
	Nx和Ny分别表示了在X方向和Y方向上划分的点数，因而网格数需要加1
	"""
	dx = (Xmax - Xmin)/(Nx+1)
	dy = (Ymax - Ymin)/(Ny+1)
	D, H = Generate_Diagonal_Matrix(Nx, Ny, dx, dy)
	A_shape = (Nx*Ny, Nx*Ny)
	b = np.zeros(Nx*Ny)
	Border_Up, Border_Down, Border_Left, Border_Right = Border_Generator(Nx, Ny)
	for i in range(Nx):
		for j in range(Ny):
			if 0 == i:
				b[Ny*i+j] = b[Ny*i+j] + (-dy**2)*Border_Left[j+1]
			if (Nx-1) == i:
				b[Ny*i+j] = b[Ny*i+j] + (-dy**2)*Border_Right[j+1]
			if 0 == j:
				b[Ny*i+j] = b[Ny*i+j] + (-dx**2)*Border_Up[i+1]
			if (Ny-1) == j:
				b[Ny*i+j] = b[Ny*i+j] + (-dx**2)*Border_Down[i+1]
	b = np.mat(b).T
	return A_shape, b, D, H

def Generate_Diagonal_Matrix(Nx, Ny, dx, dy):
	Block_Line_String = str(dx**2)+","+str((-2)*(dx**2+dy**2))+","+str(dx**2)+",0"
	#-3是因为在上一行已经产生了4个元素，而第1个元素是为了生成方便加上的，因此在上一行产生了3个元素
	#这里也可以看出，D的维度是由Ny决定的
	Block_Line_String = Block_Line_String + ",0"*(Ny-3)
	index = Block_Line_String.find(",")+1
	Diagonal_Matrix_String = Block_Line_String[index:]
	position = Block_Line_String.rfind(",")
	for i in range(Block_Line_String.count(",")-1):
		Diagonal_Matrix_String = Diagonal_Matrix_String+";"+"0,"*i+Block_Line_String[0:position]
		position = Block_Line_String.rfind(",",0,position)
	D = np.matrix(Diagonal_Matrix_String)
	Kernal_Dim = D.shape[0]
	H = np.matrix(np.eye(Kernal_Dim))*(dy**2)
	return D, H

def Border_Generator(Nx, Ny, Border_Type="array"):
	Border_Up = np.linspace(0,0,Nx+2)
	Border_Down = np.linspace(5,5,Nx+2)
	Border_Left = np.linspace(0,5,Ny+2)+np.cos(np.linspace(0,5,Ny+2))
	Border_Right = np.linspace(0,5,Ny+2)+np.sin(np.linspace(0,5,Ny+2))
	if Border_Type == "array":
		return Border_Up, Border_Down, Border_Left, Border_Right
	if Border_Type == "matrix":
		return np.mat(Border_Up[1:-1]), np.mat(Border_Down[1:-1]), np.mat(Border_Left).T, np.mat(Border_Right).T

